Crystalline boson phases in harmonic traps: 
Beyond the Gross-Pitaevskii mean field 
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Strongly repelling bosons in two-dimensional harmonic traps are described through breaking of 
rotational symmetry at the Hartree-Fock level and subsequent symmetry restoration via projection 
techniques, thus incorporating correlations beyond the Gross-Pitaevskii (GP) solution. The bosons 
localize and form polygonal-ring-like crystalline patterns, both for a repulsive contact potential and 
a Coulomb interaction, as revealed via conditional-probability-distribution analysis. For neutral 
bosons, the total energy of the crystalline phase saturates in contrast to the GP solution, and its 
spatial extent becomes smaller than that of the GP condensate. For charged bosons, the total energy 
and dimensions approach the values of classical point-like charges in their equilibrium configuration. 
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PACS numbers: 03.75.Hh, 03.75.Nt 

Bose-Einstein condensates (BEC's) in harmonic traps 
1] are normally associated with weakly interacting neu- 
tral atoms, and their physics is described adequately by 
the Gross-Pitaevskii (GP) mean-field theory Lately, 
however, experimental advances in controlling the inter- 
action strength [3J, |J, |5|, |(j permit the production of novel 
bosonic states in the regime of strong interparticle repul- 
sions. Theoretical efforts motivated by this capability 
include studies of the Bose-Hubbard model 0, and 
investigations about the "fermionization" limit of an one- 
dimensional (ID) gas of trapped impenetrable bosons 
[HI [13. often referred to as the Tonks- Girardeau (TG) 
regime [9l fl2) . In this Letter, we address the still open 
problem of strongly repelling (impenetrable) bosons in 
higher dimensions, and in particular in two dimensions 
(2D). 

We describe the strongly repelling bosons through 
symmetry breaking at the Hartree-Fock (HF) mean-field 
level followed by post-Hartree-Fock symmetry restora- 
tion, thus incorporating correlations beyond the GP solu- 
tion. This two-step method, which has not been applied 
yet to the bosonic many-body problem, has been shown 
to successfully describe strongly correlated electrons in 
2D semiconductor quantum dots [l3|. We focus here on 
results for 2D interacting bosons in a harmonic trap, with 
the extension to 3D systems being straightforward. 

To illustrate our method, we consider systems with a 
few bosons. The method describes the transition from a 
BEC state to a crystalline phase, in which the trapped 
localized bosons form crystalline patterns. In 2D, these 
patterns are ring-like, both for a repulsive contact and 
a Coulomb interaction. At the mean-field level, these 
crystallites are static and are portrayed directly in the 
single-particle densities. After restoration of symmetry, 
the single-particle densities are rotationally symmetric, 
and thus the crystalline symmetry becomes "hidden" ; 
however, it can be revealed in the conditional proba- 
bility distribution (CPD, anisotropic pair correlation), 
P(r, ro), which expresses the probability of finding a par- 



ticle at r given that the "observer" (i.e., reference point) 
is riding on another particle at ro |l4| . 

Mean-field symmetry breaking for bosonic systems 
has been discussed earlier in the context of two- 
component condensates, where each species is asso- 
ciated with a different space orbital hj. We con- 
sider here one species of bosons, but allow each par- 
ticle to occupy a different space orbital </>i(rj). The 
permanent \$n) — Pem[0i(ri), 0jv(rjy)] serves as 
the many-body wave function of the unrestricted Bose- 
Hartree-Fock (UBHF) approximation. This wave func- 
tion reduces to the Gross-Pitaevskii form with the 
restriction that all bosons occupy the same orbital 
00 (r), i-e., |$^ p ) = IIiIi < A)(r l ), and (j> (r) is deter- 
mined self-consistentlv at the restricted Bose-Hartree- 
Fock (RBHF) level via the equation 17] [HoOi) + 
(N - l)j dv 2 ^{v 2 )V(v 1 ,v 2 ) ( j )Q ^2)]M^) = eo0o(ri). 
Here V(yi,y 2 ) is the two-body repulsive interaction, 
which can be either a long-range Coulomb force, Vc = 
Z 2 e 2 /(«|ri — r2 1), for charged bosons or a contact poten- 
tial, Vs = gS(ri — r%), for neutral bosons. The single- 
particle hamiltonian is given by Ho(r) = — ?i 2 V 2 /(2m) + 
muj^r 2 /2, where ujq characterizes the harmonic confine- 
ment. 

First step: Symmetry breaking. Going beyond the 
GP approach to the unrestricted Hartree-Fock level (i.e., 
using the permanent |$jv)) results in a set of UBHF 
equations with a higher complexity than that encoun- 
tered in electronic structure problems [13(a)]. Conse- 
quently, we simplify the UBHF problem by consider- 
ing explicit analytic expressions for the space orbitals 
<j)i{Yi). In particular, since the bosons must avoid oc- 
cupying the same position in space in order to mini- 
mize their mutual repulsion, we take all the orbitals 
to be of the form of displaced Gaussians 0]> namely, 
0i(rj) = Tr -1 / 2 ^ -1 exp[— (rj - a l ) 2 /(2a 2 )]. The positions 
a; describe the vertices of concentric regular polygons, 
with both the width a and the radius a — \an | of the regu- 
lar polygons determined variationally through minimiza- 
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FIG. 1: Total energies as a function (a) of Rg and (b) of Rw 
for various approximation levels, calculated for N = 6 har- 
monically confined 2D bosons in the (1,5) lowest-energy con- 
figuration. Notation: RBHF/G - Restricted Bose-Hartree- 
Fock (RBHF) energy, E%bhf i with tne common orbital 
4>o{r) approximated by a Gaussian centered at the trap ori- 
gin; GP - the Gross-Pitaevskii energy; PRJ - the energy of 
the symmetry-restored state obtained via projection of the 
(unrestricted) UBHF state. Energies in units of hioo- 

tion of the total energy E'ubhf = ($n\H |<J>at) /{$n\$n}, 
where H = £\ =1 Ho(i*i) + Y,i<j V ( T i> T j) is the many- 
body hamiltonian. 

With the above choice of localized orbitals, the un- 
restricted permanent |$jv) breaks the continuous rota- 
tional symmetry. However, for both the cases of a con- 
tact potential and a Coulomb force, the resulting energy 
gain becomes substantial for stronger repulsion. Con- 
trolling this energy gain (the strength of correlations) is 
the ratio R between the strength of the repulsive poten- 
tial and the zero-point kinetic energy. Specifically, for 
a 2D trap, one has R$ = gm/{2nTi 2 ) for a contact po- 
tential and Rw — Z 2 e 2 / {TlujqIq) for a Coulomb force, 
with Iq = yh] (mwrj) being the characteristic harmonic- 
oscillator length. (The subscript W in the case of a 
Coulomb force stands for "Wigner", since the Coulomb 
crystallites in harmonic traps are finite-size analogs of 
the bulk Wigner crystal flif.) 

In Fig. 1, we display as a function of the parameters 
Rs (a) and Rw (b), respectively, the total energies for 
N = 6 bosons calculated at several levels of approxima- 
tion. In both cases the lowest UBHF energies correspond 
to a (1,5) crystalline configuration, namely one boson is 
at the center and the rest form a regular pentagon of ra- 



dius a. Observe that the GP total energies are slightly 
lower than the -Erbhf ones ; however, both exhibit an 
unphysical behavior since they diverge as Rs —> oo. This 
behavior contrasts sharply with that of the unrestricted 
Hartree-Fock energies, -Etjbhf, and those of the projected 
(PRJ) states (see below), which saturate as Rs — ► oo; in 
fact, a value close to saturation is achieved already for 
Rs {Rw) ~ 10- We have checked that for all cases with 
N = 2 — 7, the total energies exhibit a similar behav- 
ior. For a repulsive contact potential, the saturation of 
the UBHF energies is associated with the ability of the 
trapped bosons (independent of N) to minimize their mu- 
tual repulsion by occupying different positions in space, 
and this is one of our central results. For N = 2, the two 
bosons localize at a distance 2a apart to form an antipo- 
dal dimer. For N < 5 the preferred UBHF crystalline 
arrangement is a single ring with no boson at the center 
[usually denoted as (0, N)]. N — 6 is the first case having 
one boson at the center [designated as (1,N — 1)], and 
the (0,6) arrangement is a higher energy isomer. 

The saturation found here for 2D trapped bosons inter- 
acting through strong repelling contact potentials is an 
illustration of the "fermionization" analogies that appear 
in strongly correlated systems in all three dimensionali- 
ties. Indeed such energy saturation has been shown for 
the TG ID gas 0, 0], and has also been discussed for 
certain 3D systems (Le^ three trapped bosons [i(J and 
an infinite boson gas |2jJ). Saturation of the energy and 
the length of the trapped atom cloud (and thus of the in- 
terparticle distance) has been measured recently for the 
ID TG gas (see in particular Fig. 3 and Fig. 4 in Ref. 
and compare to the similar trends predicted here for the 
2D case in Fig. 1 and Fig. 2). 

For the Coulomb potential [see Fig. 1(b)], the dis- 
played total energies have been referenced to the clas- 
sical energy Eq [22| (plus the zero-point energy) of six 
trapped point charges in their (1,5) equilibrium configu- 
ration, since the total energy of a Wigner crystallite (in- 
dependently of whether it consists of bosons or fermions) 
is expected to approach Eq as Rw —* oo. We see again 
that the ^bhf 

energies (one common Gaussian orbital) 
diverge as Rw —> oo. In contrast, the unrestricted HF 
energies -Etjbhf remain finite and approach slowly Eq as 
Rw — > oo. A similar behavior is exhibited by the total 
energies for all TV = 2 — 7 cases of charged bosons. 

In Fig. 2, we display for the N = 6 bosons the radii 
of the polygonal rings a and widths a of the Gaussian 
orbitals obtained in various approximations, as a func- 
tion of Rs (a) and Rw (b) . For the contact potential, in 
the RBHF/G approximation we find that a — and the 
width (marked as RBHF/G in Fig. 2) keeps increasing 
continuously as Rs — > oo (this reflects the unsuccessful 
attempt of the common orbital to minimize the mutual 
repulsion between the bosons by spreading out as far as 
possible). In contrast, the unrestricted widths oubhf as " 
sociated with the displaced Gaussian orbitals (that cor- 




FIG. 2: Variationally determined widths (a) and ring radii 
(a) for TV = 6 harmonically confined 2D bosons as a func- 
tion of (a) Rs and (b) Rw, obtained according to the various 
approximations (as marked in the figure, see also caption of 
Fig. 1 ). The saturation values of a of the lowest-energy con- 
figuration for 2 < N < 7 are on the right in (a). Lengths 
in units of lo- For the UBHF case [displaying an (1,5) crys- 
tallite] the interparticle distance on the pentagonal shell is 
d — ((5 — 5 1 ^ 2 )/2) 1 ' 2 o ~ 1.176a, showing the same saturation 
trend as the radius a. 

respond to a lower total energy, see Fig. 1) saturate to a 
constant value. Similar behaviors are also exhibited by 
<7 rbhf an( i o'ubhf in the case of a Coulomb force [see of 
Fig. 2(b)]. 

The radii a associated with the pentagonal ring of lo- 
calized orbitals, however, exhibit a different behavior de- 
pending on whether the repulsive potential is a contact 
or a Coulomb one. Indeed, in the Coulomb case, the radii 
Q-ubhf keep increasing with Rwi approaching the equi- 
librium radius a,Q = 1.334:loR w 3 of six Ze classical point- 
charges in a harmonic trap in the (1,5) configuration |22j . 
In contrast, for a repulsive contact potential, the radii 
otjbhf saturate to a constant value « 2lg. The depen- 
dence of the saturation values of a on TV (for 3 < N < 7) 
for the lowest-energy configurations is shown on the right 
in Fig. 2(a). The different behavior of the boson positions 
in the UHBF crystallite is a natural consequence of the 
long-range character of the Coulomb potential versus the 
short-range contact potential. 

Second step: Restoration of broken symmetry. Al- 
though the optimized UBHF permanent |$jv) performs 
exceptionally well regarding the total energies of the 
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FIG. 3: (a-c): Single-particle densities for N = 6 2D har- 
monically trapped neutral bosons with a contact interaction 
and Rg — 25. (a) The single-orbital self-consistent GP case, 
(b) The symmetry broken UBHF case (static crystallite), (c) 
The projected (symmetry-restored wave function, see Ref. 
(2*3^ case (collectively fluctuating crystallite). The crystalline 
structure of the outer ring in this last case is "hidden", but 
it is revealed in the conditional probability distribution 
displayed in (d), where the observation point is denoted by a 
black dot (on the right). Lengths in units of Iq. 



trapped bosons, in particular in comparison to the resc- 
tricted wave functions (e.g., the GP anzatz), it is still 
incomplete. Indeed, due to its localized orbitals, |$j\r) 
does not preserve the circular (rotational) symmetry of 
the 2D many-body hamiltonian H. Instead, it exhibits 
a lower point-group symmetry, i.e., a C2 symmetry for 
N = 2 and a C5 one for N = 6 (see below). As a re- 
sult, |3>jv) does not have a good total angular momen- 
tum. This is resolved through a post-Hartree-Fock step 
of restoration of broken symmetries via projection tech- 
niques [13(b)], , yielding a new wave function I^^Y) 
|23j with a definite angular momentum L. Here, we fo- 
cus on the properties of the ground state, i.e., L = 0; the 
corresponding energy is i?^ R J . 

For N = 6 2D bosons, Fig. 1 shows that the E% RJ 
energies share with the UBHF ones the saturation prop- 
erty for the case of a contact-potential repulsion, as well 
as the property of converging to E^i as Rw -> 00 for 
the case of a Coulomb repulsion. In both cases, however, 
the projections bring further lowering [24| of the total 
energies compared to the UBHF ones. Thus, for strong 
interactions (large values of Rs or Rw) the restoration- 
of-broken-symmetry step yields an excellent approxima- 
tion of both the exact many-body wave function and the 
exact total energy |25j . 

The transformations of the single-particle densities 
(displayed in Fig. 3 for N = 6 neutral bosons interact- 
ing via a contact potential and Rs = 25) obtained from 
application of the successive approximations provide an 
illustration of the two-step method of symmetry breaking 
with subsequent symmetry restoration. Indeed, the GP 
single-particle density [Fig. 3(a)] is circularly symmet- 
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ric, but the UBHF one [Fig. 3(b)] explicitly exhibits a 
(1,5) crystalline configuration. After symmetry restora- 
tion [Fig. 3(c)], the circular symmetry is re-established, 
but the single-particle density is radially modulated un- 
like the GP density. In addition, the crystalline struc- 
ture in the projected wave function is now hidden; how- 
ever, it can be revealed through the use of the CPD [l^j 
[see Fig. 3(d)], which resembles the (crystalline) UBHF 
single-particle density, but with one of the humps on 
the outer ring missing (where the observer is located). 
In particular, F(r ,r ) ~ and the boson associated 
with the observer is surrounded by a "hole" similar to 
the exchange-correlation hole in electronic systems. This 
is another manifestation of the "fermionization" of the 
strongly repelling 2D bosons. However, here as in the ID 
TG case |9|,ll2|, the vanishing of -P(i"o, ro) results from the 
impenetrability of the bosons. For the GP condensate, 
the CPD is independent of r , i.e., PGp( r , r o) oc |</>o(r)| 2 , 
reflecting the absence of any space correlations. 

It is of importance to observe that the radius of the 
BEC [GP case, Fig. 3(a)] is significantly larger than the 
actual radius of the strongly-interacting crystalline phase 
[projected wave function, Fig. 3(c)]. This is because the 
extent of the crystalline phase saturates, while that of 
the GP condensate grows with no bounds as Rs — > oo. 
Such dissimilarity in size (between the condensate and 
the strongly-interacting phase) has been also predicted 
[To) for the trapped ID Tonks- Girardeau gas and indeed 
observed experimentally |(J. In addition, the 2D single- 
particle momentum distributions for neutral bosons have 
a one-hump shape with a maximum at the origin (a be- 
havior exhibited also by the trapped ID TG gas). The 
width of these momentum distributions versus Rs in- 
creases and saturates to a finite value, while that of the 
GP solution vanishes as Rs — > oo. 

In conclusion, we provided a solution to strongly re- 
pelling bosons in 2D harmonic traps using a two-step 
method of breaking of rotational symmetry at the un- 
restricted Bose-Hartree-Fock level and of subsequent 
symmetry restoration. This method yields substan- 
tially lower total energies compared to the GP solu- 
tion, through the inclusion of correlations beyond the 
single-orbital Bose-Einstein condensate. We find that 
the bosons become localized and form crystalline pat- 
terns made of concentric polygonal rings, both for a re- 
pulsive contact and a Coulomb interaction. For neutral 
bosons the total energy of the crystalline phase saturates 
with increasing strength of the repulsion, in contrast to 
the GP condensate whose energy diverges. Furthermore, 
the spatial extent saturates and becomes smaller than 
that of the GP condensate, which grows without limit. 
For charged bosons, the total energy and spatial extent 
of the crystalline phase approach the classical values of 
point-like charges in their equilibrium configuration as 
Rw — > oo. In light of the above, we trust that our pre- 
dictions will provide the impetus for experimental efforts 



to access the regime of strongly repelling bosons in two 
dimensions. To this end we anticipate that extensions of 
methodologies developed for the recent realization of the 
Tonks-Girardeau regime in ID (using a finite small num- 
ber of trapped 87 Rb and optical lattices, with a demon- 
strated wide variation of Rg from 5 to 200 |f| and from 
1 to 5 0) will prove most promising. Control of the in- 
teraction strength via the use of the Feshbach resonance 
may also be considered Q. 
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